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Abstract 

Background: Water temperature greatly influences the physiology and behaviour of teleost fish as other aquatic 
organisms. While fish are able to cope with seasonal temperature variations, thermal excursions outside their 
normal thermal range might exceed their ability to respond leading to severe diseases and death. 
Profound differences exist in thermal tolerance across fish species living in the same geographical areas, promoting 
for investigating the molecular mechanisms involved in susceptibility and resistance to low and high temperatures 
toward a better understanding of adaptation to environmental challenges. The gilthead sea bream, Spams aurata, is 
particularly sensitive to cold and the prolonged exposure to low temperatures may lead to the "winter disease", a 
metabolic disorder that significantly affects the aquaculture productions along the Northern Mediterranean coasts 
during winter-spring season. While sea bream susceptibility to low temperatures has been extensively investigated, 
the cascade of molecular events under such stressful condition is not fully elucidated. 

Results: In the present study two groups of wild sea bream were exposed for 21 days to two temperature regimes: 
16±0.3°C (control group) and 6.8±0.3°C (cold-exposed group) and DNA microarray analysis of liver transcriptome 
was carried out at different time points during cold exposure. 

A large set of genes was found to be differentially expressed upon cold-exposure with increasingly relevant effects 
being observed after three weeks at low temperature. All major known responses to cold {i.e. anti-oxidant response, 
increased mitochondrial function, membrane compositional changes) were found to be conserved in the gilthead 
sea bream, while, evidence for a key role of unfolded protein response (UPR) to endoplasmic reticulum (ER) stress, 
during short- and long-term exposure to cold is reported here for the first time. 

Conclusions: Transcriptome data suggest a scenario where oxidative stress, altered lipid metabolism, ATP depletion 
and protein denaturation converge to induce ER stress. The resulting UPR activation further promotes conditions for 
cell damage, and the inability to resolve ER stress leads to severe liver dysfunction and potentially to death. 
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Background 

Teleosts are ectotherms and water temperature greatly 
influences their physiology and behaviour. Within the 
species-specific thermal tolerance ranges, fishes are able 
to cope with natural temperature changes, such as daily 
and seasonal variations. However, close to or beyond the 
limits of thermal tolerance, physiological adaptation is 
no longer capable to cope with environmental condi- 
tions and fish health and survival are threatened [1,2]. 
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Thermal stress trigger cellular and organismal responses 
activating a cascade of events that may lead to patho- 
logical conditions and ultimately to death. There are, 
however, profound differences in thermal tolerance even 
among fish species with similar geographic distribution 
[1-3]. Investigating the molecular mechanisms under- 
lying the susceptibility to cold temperatures could pro- 
vide a better understanding of how organisms either 
adapt or fail to respond to environmental challenges. 
The gilthead sea bream, Spams aurata Linneus 1758, al- 
though a relatively robust species, is particularly sensi- 
tive to cold: when water temperature drops below 12°C, 
fish become minimally active, reducing food intake and 
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metabolism. For this species the lethal temperature is 
higher (around 5°C) compared to other teleosts with 
similar ecological traits and geographic distribution {e.g., 
European sea bass Dicentrarchus labrax, meagre Argyro- 
somus regius) [3,4]. 

In the gilthead sea bream, the persistence of low temper- 
atures for a long time may predispose to the "winter dis- 
ease", a syndrome that involves several factors such as 
generalized stress, metabolic depression, immune suppres- 
sion and opportunistic pathogen infections. Along the 
Northern Mediterranean coasts and lagoons, the outbreak 
of winter disease may lead to fish mass mortality up to 
80% during the winter season, with significant economic 
losses for semi-intensive and extensive aquaculture pro- 
duction [5-9]. 

Several aspects of the physiological response to low tem- 
peratures in the gilthead sea bream have been investigated 
both in natural populations and under laboratory condi- 
tions [3,10-13], providing evidence that metabolic depres- 
sion and oxidative stress occur in different tissues after cold 
exposure. One of the main pathological consequences of 
the limited response to cold exposure in the sea bream is 
the onset of a severe hepatic steatosis, suggesting that liver 
dysfunction might be at the basis of several winter disease 
symptoms [10]. However, the cascade of molecular events 
underlying the response to cold in this temperature- 
sensitive species is not fully elucidated. Winter disease rep- 
resents an extreme example of the complex interactions 
between organismal physiology and environmental stressors, 
which remain a central issue in aquatic animal ecology [14]. 
Until recently, biochemical and molecular studies on organ- 
ismal adaptation to environmental variations were limited to 
single or few proteins or genes. High-throughput transcrip- 
tome analysis, however, now provides an unprecedented view 
of global transcriptome response to environmental stressors. 
In particular, DNA microarray and next-generation se- 
quencing, provided detailed molecular information on 
genomic responses of teleosts to hypoxia [15,16], heat 
shock [17-20], cold shock/cold stress [16,21-28]. In the 
present study an oligo-DNA microarray representing 
19,7345. aurata unique transcripts [29] was used to inves- 
tigate the response to cold on liver, as this organ exerts 
important metabolic functions in cold adaptation pro- 
cesses [10,21]. Fish were exposed for 21 days to 6.8°C in 
order to simulate environmental conditions that frequently 
occur during winter in the Northern Mediterranean coastal 
areas. This critical temperature is near to the lethal thresh- 
old reported for this species. The experiment was carried 
out during winter season, to prevent seasonal impact that 
might influence the physiological status of the fish. DNA 
microarray analyses were carried out at four different time 
points (0 h, 6 h, 24 h and 21 day) in order to compare liver 
transcriptome profiles in cold-exposed and control individ- 
uals (16 ± 0.3°C). 



Results 

Growth and liver condition 

The effects of low temperature exposure on growth and 
liver condition of sea bream are summarized in Table 1. 
The body weight (BW) was only slightly reduced in cold- 
exposed groups at the end of the experiment (21d), with no 
differences observed at starting time (0 h) with controls. 

Conversely, liver weight and hepatosomatic index (HSI) 
significantly increased in cold-exposed groups (6.8 ± 0.3°C) 
at 21d and in comparison with the control groups (16 ± 
0.3°C). In most of fish exposed to cold, the liver appeared 
yellowish and friable and affected by steatosis. No mortal- 
ity occurred during the trial. 

Microarray analysis 

A Principal Component Analysis (PCA) of microarray 
data (Figure 1) showed that 72.4% of the observed tran- 
scriptional changes were explained by the first two com- 
ponents and 6.2% by the third component. Since along 
the first component samples are ordered according to 
time of exposure, it seems evident already from such an 
exploratory analysis that long-term exposure to cold 
(21 days) caused the largest transcriptional variation be- 
tween cold-exposed and control groups. 

To evaluate the effects of low temperature, pair-wise 
comparisons between cold-exposed and controls at each 
time point using a two-class SAM test (FDR = 5%; FC > 2; 
see Methods) were carried out, showing large sets of over- 
and under-expressed genes in cold-exposed animals, with 
increasingly relevant effects being observed at later time 
points (see Table 2). A total of 198 (109 up- and 89 under- 
expressed), 1,419 (709 over- and 710 under-expressed), 
3,124 (1,548 over- and 1,576 under-expressed) and 4,194 
(1,835 over- and 2,359 under-expressed) differentially 
expressed genes (DEGs) were obtained at 0 h, 6 h, 24 h 
and 2 Id, respectively (Additional file 1). A substantial over- 
lap was observed between DEGs at different time points 
(see Figure 2), although just 72 transcripts were signifi- 
cantly differentially expressed along the entire experiment. 
The response to cold was strong not only in terms of total 

Table 1 Mean value (±SD) of growth parameters and liver 
condition in control (CTRL) and cold exposed groups (COLD) 
at the start (0 h) and at the end (21 d) of experiment 

Time Oh 21 d 

CTRL COLD CTRL COLD 

BW(g) 121.3±17.3 a 121.2±14.7 a 118.6±17.5 a 119.7±20.3 a 

TL (cm) 21.2±1.0 a 21.3±0.9 a 21.5 ±1.1 3 21.3±1.0 a 

LW(g) 1 .87 ± 0.72 a 2.07 ± 0.63 a 1.10±0.45 b 2.48±0.85 c 

HSI 1 .52 ± 0.48 a 1.7±0.44 b 0.92 ± 0.33 c 2.04 ± 1.5 1 d 

Steatosis (%) 25 a 35 a 5 a 100 b 

Different letters represent significant differences between groups and 
sampling time at p<0.05. 
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Figure 1 PCA of S. aurata liver gene expression profiles in controls (circle) and cold groups (square) at the different sampling points 



number of DEGs, but also in the degree of transcriptional 
changes. Thirty DEGs at 0 h (21 over- and 9 under- 
expressed), 202 at 6 h (141 over- and 61 under-expressed), 
485 at 24 h (240 over- and 245 under-expressed), 660 at 
21d (318 over- and 342 under-expressed) showed large 
FCs (log 2 FC>2) (see Table 2). The time component in 

Table 2 Absolute number of differentially expressed 
genes (DEGs) at each time-points after temperature drop 
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Figure 2 Venn diagram showing the number of common 


differentially expressed genes at the four sampling time-points, 


obtained by SAM. 
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transcriptome regulation appears to be largely due to the 
effect of cold exposure, as SAM tests performed compar- 
ing only control samples across different time points 
showed no differentially expressed genes (FDR = 5%). Such 
evidence confirms that the liver transcriptome in control 
animals is relatively unaffected by the duration of the ex- 
periment per se, despite the fact that controls were also 
fasted to increase comparability with cold-exposed ani- 
mals. A complete list of DEGs is available in Additional 
file 1. 

A useful tool for reducing the complexity of large sets of 
DEGs is functional annotation and enrichment analysis, 
which highlights the most significant biological processes 
across DEGs (Additional file 2). The most significant 
process at 0 h was G 0:0043433 ~ negative regulation of 
transcription factor activity, with three DEGs: heme oxy- 
genase (decycling) 1 (HMOX1), CAAT /enhancer binding 
protein gamma (CEBPG), and DNA-damage-inducible 
transcript (DDIT3), more commonly reported as CHOP. 
All three genes have been reported to be involved in anti- 
oxidant response, although DDIT3/CHOP is best known 
as one of the main effectors of the Unfolded Protein Re- 
sponse (UPR) [30]. Highly relevant to response to oxida- 
tive stress, is also Kelch-like ECU- associated protein 1 
(KEAP1, Additional file 1), a cytosolic inhibitor of nuclear 
factor erythroid 2-like 2 (NFE2L2, better known as NRF2). 
Changes in the cellular redox state induce NFE2L2/NRF2 
dissociation from KEAP1, its nuclear import and the 
transcriptional activation of several antioxidant genes, 
including HMOX1, and of its own inhibitor, KEAP1. 
Up-regulation of KEAP1 and HMOX1, as observed at 
0 h as well as 6 h, 24 h and 2 Id, suggests an immediate 
and sustained activation of NFE2L2/NRF2, the main 
transcription factor involved in the cellular response to 
oxidative stress. 

At 6 h, the most significant enriched Biological Process 
was GO:0044255 ~ cellular lipid metabolic process repre- 
sented by 31 DEGs (Additional file 2). Among them, it is 
listed Glycerol-3-phosphate acyltransferase mitochondrial 
precursor (GPAT-1), that represents the first step of the 
metabolic pathway for the synthesis of glycerolipids [31]. 
Its transcript was significantly over-expressed at 6 h, but 
also at 24 h and 2 Id. Similar evidence (over-expression at 
6 h, 24 h, and 2 Id) was observed for acyl-CoA synthetase 
long-chain family member 4 (ACSL4), a gene that was 
reported to be over-expressed in nonalcoholic fatty 
liver (NAFL) [32]. The opposite pattern (constant down- 
regulation at 6 h-21d) was observed for other genes in the 
list, acyl-CoA dehydrogenase C-2 to C-3 short chain 
(ACADS), acyl-CoA dehydrogenase long chain (ACADL), 
and enoyl-coenzyme A hydratase, all playing important 
roles in lipid beta-oxidation. On the same line of evidence 
is the observed under-expression at 6 h-21d of hydroxy ste- 
roid (17-beta) dehydrogenase 4 (HSD17B4), which is 



involved in the peroxisomal pathway of lipid catabolism. 
Another gene belonging to the cellular lipid metabolic 
process is, noteworthy, peroxisome proliferator-activated 
receptor delta (PPARD), which is significantly under- 
expressed at 6 h and also at 2 Id. Activation of PPARD was 
recently reported to attenuate hepatic fat deposition [33] . 

Finally, a sea bream transcript encoding fatty acid binding 
protein 6 (FABP6) was significantly under-expressed at 6 h, 
as well as at all other time points (Additional file 1). At least 
three other transcripts encoding FABPs (FABP1, FABP2; 
FABP7) were found to be under- or over-expressed, espe- 
cially at 21d (see Figure 3). The most relevant for fatty acid 
transport within hepatic cells is FABP1, which is the main 
FABP expressed in the liver [34], as confirmed by absolute 
fluorescence values in this study (Figure 3). Multiple lines 
of evidence support the hypothesis that FABP1 acts as a 
long-chain fatty acid transporter, specifically targeting its li- 
gands to beta-oxidation pathways [34]. In the liver of cold- 
exposed sea breams FABP1 drops dramatically after 21 days 
(Figure 3). 

In consideration of the extremely large number of differ- 
entially expressed genes at 24 h and 2 Id, enrichment ana- 
lysis was performed separately for over- and under- 
expressed DEGs (Additional file 2), although the results 
were similar when considering all DEGs together (data 
not shown). At 24 h, the most significantly enriched 
biological process (BP) for over- expressed genes was 
GO:0006396 ~ RNA processing with 47 DEGs, the major- 
ity of which represented proteins involved in RNA splicing 
(Additional file 2). In fact, the second most significant BP 
was GO:0008380 ~ RNA splicing with 34 transcripts. 
Other noteworthy BPs were GO:0046907 ~ intracellular 
transport and GO:0008219 ~ cell death. In the latter BP it 
is included Growth arrest and DNA damage-inducible 
beta (GADD45B) gene, which plays a key role in the re- 
sponse to cellular stress and it has been recently reported 
to be transcriptionally activated by oxidative stress in hep- 
atic cells [35]. While GADD45B has an essential pro- 
survival role, other over-expressed genes have opposite 
actions. For instance, BCL2-like 14 (BCL2L14) belongs to 
the pro-apoptotic branch of the Bcl2 family, cell death- 
inducing DFFA-like effector c (CIDEC), phosphatase and 
tensin homolog (PTEN), tumor protein p53 binding protein 
2 (TP53BP2), all favor apoptosis at different levels, while 
caspase 6 (CASP6) is one of executioner enzymes for pro- 
grammed cell death. In the same list is also the transcript 
encoding CHOP, which was found to be over-expressed as 
early as 0 h and is among the genes that are differentially 
expressed at all time points (Additional file 1). In fact, per- 
sistent over-expression of CHOP is considered a marker of 
chronic and/or excessive activation of UPR, and CHOP has 
been shown to promote cell death [36]. Particularly, the 
UPR branch leading to CHOP transcription has been shown 
to be involved in diverse liver diseases (e.g fatty liver disease 
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Figure 3 Average of gene expression values of FABP1, FABP2, FABP6 and FABP7 at each sampling time (0 h = 0 hours; 6 h 6 hours; 
24 h 24 hours; 21d = 21 days) in control and cold-exposed groups. 



(FLD) [36,37]). The most significant BPs for under-expressed 
genes at 24 h concerned lipid metabolism (GO:0044255 ~ 
cellular lipid metabolic process, GO:0032787 ~ monocarbox- 
ylic acid metabolic process, GO:0006629 ~ lipid meta- 
bolic process, GO:0008610 ~ lipid biosynthetic process). 
Several under-expressed transcripts encoded the same 
proteins already described for 6 h (ECHDC2, ACSS2, 
HSD17B4), and additional ones involved in lipid beta- 
or alpha-oxidation, (e.g. carnitine O-acetyltransf erase, 
enoyl Co A hydratase 1, peroxisomal hydroxy acid oxi- 
dase 2), or in the synthesis of phospholipids (glycer one- 
phosphate O-acyltransferase, phosphatidylethanolamine 
N-methyltransf erase) . 

Gene expression profiles were even more divergent be- 
tween controls and cold-exposed fish after 21 days at con- 
stant low temperature. While several of BPs enriched at 
earlier time points remain significant (e.g. lipid meta- 
bolism, RNA processing), the most significantly over- 
represented pathways/processes involve protein metabol- 
ism, in particular protein catabolism and folding. Such 
evidence lends further support to the hypothesis of UPR 
activation in sea bream hepatic cells under cold stress. In 
fact, Gene set enrichment analysis (GSEA; see methods) 
using the reactome "Unfolded Protein Response (UPR)" 



gene set was marginally significant at 24 h (FDR < 25%) 
and highly significant at 21d (adjusted p < 0.01) (Additional 
file 3). In addition to the mentioned CHOP, GSEA identi- 
fied as significantly over-expressed in cold-exposed fish 
two other key UPR effectors, X-box binding protein 1 
(XBP1) and activating transcription factor 4 (ATF4). It has 
been recently shown that prolonged activation of ATF4 
and CHOP leads to cell death through a specific mechan- 
ism [38]. While UPR normally decreases protein synthesis, 
joint over-expression of ATF4 and CHOP activates protein 
synthesis-related genes, in addition to the UPR trans- 
criptional program. Increased protein synthesis, with 
consequent ATP depletion, coupled with enhanced oxida- 
tive stress, leads to cell death [38]. Remarkably, in cold- 
exposed animals at 21d nine aminoacyl-tRNA synthetase 
(ARS) genes (YARS, WARS, TARS, AARS, VARS, SARS, 
GARS, IARS, CARS) are over-expressed, largely overlap- 
ping with those reported to be CHOP/ATF4 common tar- 
gets by Han et al. [38]. Likewise, expression of six genes 
encoding eukaryotic translation initiation factors (EIFs) 
was significantly increased (SAM pairwise test, Additional 
file 1), including EIF2S2, EIF4G2, and EIF5, which were 
found to be activated upon CHOP/ATF4 over-expression 
in human cells. 
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GSEA also confirmed the main lines of evidence from 
the analysis of SAM-derived DEGs. When sets compris- 
ing genes included in reactomes "Triglyceride Biosyn- 
thesis" and "Mitochondrial Fatty Acid Beta Oxidation" 
were employed, GSEA revealed highly significant enrich- 
ment at 2 Id (Additional file 3). Similar evidence was ob- 
tained when analysing the set of putative target genes of 
the transcription factor CHOP (TFT: V$CHOP_01), 
while the analogous analysis on genes putatively regu- 
lated by the NRF2 (TFT: V$NRF2_01), the key regulator 
of oxidative stress response, GSEA was highly significant 
already at 6 h. With regard to proteins involved in pro- 
tection against reactive oxygen species (ROS), a gene-by- 
gene analysis revealed that classical anti-oxidant en- 
zymes such as superoxide dismutase (SOD1 and SOD2) 
are not modified at the transcriptional level, while cata- 
lase mRNA is significantly under-expressed (Additional 
file 1). On the other hand, other enzymes with similar 
function, such as peroxiredoxin 5 (PRDX5) and glutathi- 
one peroxidase 4 (GPX4) are over-expressed at 21d. Both 
PRDX5 and GPX4 have a major role in protecting mito- 
chondria from ROS damage. 

Metabolic depression and reduced oxygen consumption 
have been reported in gilthead sea breams exposed to 
temperatures lower than 12°C [39] and apparently no 
temperature compensation of metabolic rate occurs in this 
species. In eurythermal fish such compensation is nor- 
mally obtained through enhanced aerobic respiration and 
mitochondrial biogenesis. To explore whether transcripts 
involved in these processes are differentially expressed in 
cold-exposed sea breams, GSEA was carried out using two 
reactome gene sets, "TCA cycle and respiratory electron 
transport" and "Mitochondrion organization and biogen- 
esis". Genes included in the latter process were marginally 
enriched at all time points and significantly at 24 h, while 
transcripts related to aerobic respiration and tricarboxylic 
acid (TCA) cycle, were differentially expressed at 6 h-21d 
(Additional file 3). Similarly related to mitochondrial bio- 
genesis are prohibitins (PHB1 and PHB2), which have 
been shown to play a key role in mitochondrial biology 
[40]. PHB1 was found to be significantly over-expressed 
under pairwise SAM test in cold-exposed animals at 24 h 
and 21d. 

Another well-known response to cold temperature in 
fish is the over-expression of a- and ^-tubulins (TUBA 
and TUBB) and A-9 desaturase (SCD), which is involved 
in homeoviscous adaptation (HVA) [21,28]. In fact, 
TUBA6 and TUBA1 are over-expressed at all time 
points (0 h-21d), TUBB5 and TUBB2B at 0 h-21d, while 
SCD is over-expressed at 24 h and 2 Id, with a FC > 120 
at the latter time point. Likewise, mRNA encoding for 
high mobility group box 1 (HMGB1; [21]) was found to 
be significantly over-expressed at 24 h and 2 Id (Additional 
file 1). 



Validation of microarray results using Real-time RT-PCR 

Target genes for qRT-PCR analysis were selected consid- 
ering transcripts with a FC between 0.2 and 13 and rep- 
resented by two independent probes in the normalized 
data set. According to these criteria, 10 genes were se- 
lected (Additional file 4) and tested on 2 time points. 

For all target genes tested, the direction of change in 
expression was concordant between qRT-PCR data 
and microarray results. Overall, a statistically significant 
correlation was obtained comparing the expression 
levels of each gene across all biological replicates. Five 
genes showed highly significant correlation coefficients 
(Spearman rho > 0.8) for both probes (p-value < 0.01) 
with qPCR data (Additional file 4). One gene (Glycerol- 
3-phosfateacyl transferase; GPAT) exhibited a high cor- 
relation (Spearman rho > 0.8, p-value < 0.01) for one 
probe and a significant correlation (0.7 < rho >0.8 with 
p-value <0.01) for the other one. One gene (SOD) had a 
significant correlation for one probe (0.5 < rho < 0.8 
with p-value < 0.01) and a positive correlation for the 
other one (rho = 0.5, p-value > 0.01). The remaining 
three genes (26S proteasome complex subunit, malate 
dehydrogenase, ACP) showed a not significant, albeit 
positive correlation for both probes. 

Discussion 

Gene expression analyses of gilthead sea bream exposed to 
low temperature suggested the conservation of the main 
mechanisms/processes contributing to setting the thermal 
limits in ectotherms [41]. These include a broad activation 
of genes involved in RNA processing, especially at 24 h 
and 2 Id. RNA chaperones and, in general, modified RNA 
processing are thought to be required to compensate for 
increased stability of RNA secondary structures at low 
temperatures. Likewise, over-expression of SCD leads to 
increasing polyunsaturated fatty acids in lipid by-layers, 
which decrease membrane fluidity contributing to HVA. 
Up-regulation of HMGB1 was suggested to be involved in 
compensating the effects of cold on DNA [21], as HMGB1 
seems to increase "openness" of DNA to transcription on 
a genome-wide level. An active transcriptional program 
for mitochondrial biogenesis was identified as well. Signifi- 
cant alterations in the expression of mRNAs encoding pro- 
teins involved in aerobic respiration were found, although, 
as already observed in other teleosts, the net effect of such 
changes was not clear. 

The investigation of acute stress and mortality at ther- 
mal extremes might not reflect the molecular and organis- 
mal mechanisms underlying the capacity to respond to 
less extreme temperatures, which might represent the rele- 
vant thermal boundary for a certain species in nature [41]. 
The experiment presented here builds upon existing 
knowledge at different levels (biochemical, histological, 
proteomic) on a well-studied syndrome, the winter disease, 
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which represents a case study for the limited ability to re- 
spond to low temperatures. Earlier work already identified 
tissue- and organ-specific alterations as well as organismal- 
level effects of low temperatures, suggesting that metabolic 
unbalances in the liver are one of the central events for 
cold-induced disease and lethality in the gilthead sea 
bream [13]. 

Here, whole-transcriptome evidence is provided on 
acute and sub-acute reaction in such a key organ, showing 
a massive response especially after three weeks at low 
temperature. It should be noted that such a response does 
not represent the broad alterations of expression profile 
generally observed in pre-mortal conditions, as no mortal- 
ity occurred and all fish gradually recovered when water 
temperature was backed to 16°C [42]. 

Gene expression analyses highlighted the activation of 
the anti-oxidant response already at 0 h, and was sus- 
tained until 2 Id (Additional file 2 and Additional file 3). 

Recent studies in ectotherms demonstrated increased 
oxidative stress as temperature moves away from ther- 
mal optimum [41,43]. Noteworthy, the measurement of 
oxidative stress and the proteomic analyses of sea bream 
liver exposed to low temperature, revealed a significant 
increase of lipid peroxidation and nitric oxide in cold- 
exposed animals. Expression of anti-oxidant proteins 
was less clear as glutathione S-transferase and catalase 
were significantly under-expressed, in agreement with 
transcriptomic results presented here, while other en- 
zymes potentially involved in oxidative stress response 
were over-expressed [10]. A vast literature in human 
medicine suggests a relevant role for oxidative stress in 
liver disease, especially in hepatic steatosis. Such a 
knowledge is much pertinent as one of the major histo- 
logical findings in winter disease is increased lipid con- 
tent in sea bream liver [11,44]. More specifically, it has 
been reported that human hepatocytes exposed to H 2 0 2 
show decreased expression of PPARA and of its target 
genes [45] . PPARA positively regulates fatty acid oxidation, 
therefore down-regulation of PPARA and beta-oxidation- 
related transcripts as observed in human hepatocytes as 
well as in cold-exposed sea bream liver, might reduce lipid 
beta-oxidation with consequent lipid accumulation, lead- 
ing to hepatic steatosis [46,47]. This hypothesis is consist- 
ent with the liver condition of fish observed in this study. 

All the mentioned reactions {i.e. anti-oxidant response, 
increased mitochondrial function, stability-compensation 
at the RNA and DNA level, membrane compositional 
changes), however, have been already reported in cold- 
adapted teleosts or in eurythermal fish with good toler- 
ance for low water temperatures. What is different in 
the gilthead sea bream that makes it unfit for cold? The 
most interesting finding in the present study is the evi- 
dence for a key role of UPR in the acute and sub-acute 
response to cold. To our knowledge, this is the first time 



that UPR is observed after exposure to low temperatures. 
The UPR or endoplasmic reticulum (ER) stress response is 
a complex mechanism that is conserved across multicellu- 
lar organisms and is aimed at preventing ER overload with 
misfolded or damaged proteins [30]. It is particularly im- 
portant for cells that have high levels of secreted protein 
synthesis. This is the case of hepatocytes, which actively 
secrete several serum proteins {e.g. albumin, coagulation 
enzymes, complement, components), and in particular, li- 
poproteins to re-circulate lipids. UPR signaling consists of 
three main branches with diverse modes of activation {e.g. 
differential splicing, protein phosphorylation, protein mat- 
uration), which converge on a broad transcriptional pro- 
gram to counteract ER stress. Influx of newly synthesized 
proteins into the ER is decreased, ER-specific chaperones 
and proteins favoring disulfide bond formation {e.g. EROl- 
like beta, protein disulfide isomerases) are over-expressed, 
irreversibly damaged peptides are degraded by up-regulated 
ERAD components. UPR in hepatocytes has been exten- 
sively studied because of its potential relationship with liver 
disease. ER stress response represents a protective, pro- 
survival response, which is aimed at restoring cell homeo- 
stasis. UPR, however, is a double edge sword. When it is 
unable to relieve ER overload, chronic UPR activation leads 
to cell death [36]. Which are the potential triggers for UPR 
in the case of cold-exposed gilthead sea bream? Low 
temperature is known to affect protein folding. Protein de- 
naturation at 6°C might not be severe enough to comprom- 
ise protein function directly, for instance altering enzymatic 
activity, but it likely represents a stress signal for UPR. Oxi- 
dative stress might also contribute to UPR activation as free 
radicals oxidize cell proteins, and oxidized proteins are un- 
likely to pass ER quality control. Oxidative protein damage 
might be paradoxically amplified by some UPR effectors. 
For instance, EROl-like proteins (EROIL alpha and beta) 
generate ROS in the ER as part of the disulfide bond for- 
mation. Remarkably, EROlL-beta was significantly over- 
expressed at 2 Id in cold-exposed fish. Furthermore, protein 
synthesis and folding require energy. In cold-exposed sea 
breams ATP might be depleted because of lack of meta- 
bolic rate compensation [39]. Finally, prolonged fasting 
causes depletion of glycogen [12] and in the long term, re- 
duced energy production from glucose could not be com- 
pensated by fatty acid oxidation, because such metabolic 
pathway is down-regulated (PPARA and its target genes). 
As already mentioned, anti-oxidant response represses 
beta-oxidation, and UPR appears to act synergistically in 
this direction through CHOP [48]. Repression of lipid beta- 
oxidation likely contributes also to accumulation of triglyc- 
erides in sea bream hepatocytes. During cold-induced 
fasting perivisceral lipid stores are mobilized and fatty 
acids are taken up in the liver [12,49, this study]. As their 
catabolic pathway is repressed, fatty acids follow the ana- 
bolic route, which is up-regulated at the transcriptional 
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level. The observed up-regulation of PPARG and genes in- 
volved in triglyceride biosynthesis (Additional file 1 and 
GSEA results in Additional file 3) might be a response to 
increased levels of circulating lipids, although ER stress 
might contribute as well The relationships between hep- 
atic lipid biosynthesis and UPR is still controversial, al- 
though unresolved ER stress in UPR-effectors null mice is 
associated with over-expression of genes encoding key en- 
zymes for triglyceride biosynthesis, coupled with significant 
up-regulation of CHOP and ATF4 [50]. As mentioned 
earlier, prolonged over-expression of these two transcrip- 
tion factors ultimately leads to cell death through up- 
regulation of genes involved in protein synthesis, a pattern 
observed here after three weeks of cold-exposure. 

Conclusions 

Low water temperature apparently generates a perfect 
storm in the gilthead sea bream, where oxidative stress, al- 
tered lipid metabolism, ATP depletion, and protein de- 
naturation converge to induce ER stress, UPR mechanisms 
further promote conditions for cell damage, and the inabil- 
ity to resolve ER stress leads to liver dysfunction and even- 
tually to death. 

Under such a scenario, which might be the factors that 
makes the gilthead sea bream unable to cope with cold, 
compared to related fish species? Lack of temperature 
compensation for aerobic metabolism with consequent 
ATP depletion has been already proposed [39]. Limited 
capacity to prevent ROS -induced damage in this species 
might also be hypothesized. Catalase and SODs are key 
scavenger enzymes, as clearly demonstrated in several 
studies, yet SODs were unchanged and catalase was down- 
regulated at the mRNA (this study) and protein level [10]. 
Finally, UPR might be less effective in the gilthead sea 
bream, as possibly suggested by similarity in gene expres- 
sion patterns with UPR-null mice. ER stress response is 
highly conserved across the animal kingdom, but it also 
shows significant variation between natural populations in 
a single species [51]. It is therefore likely that even larger 
variation occurs between different species. 

Methods 

Animals and experimental conditions 

The experiment was performed at the Aquaculture Ex- 
perimental Centre of Veneto Agricoltura located in Valle 
Bonello (Rovigo, Italy) from January to March 2008. The 
experiment was carried out following the standards laid 
down in the Directive 98/58/EC concerning the protec- 
tion of animals kept for farming purposes, in accordance 
with the Commission Recommendation 2007/526/EC on 
guidelines for the accommodation and care of animals 
used for experimental and other scientific purposes. 

Two groups of wild sea breams (121 g body weight; 
21 cm total length) caught in Valle Bonello were kept in 



duplicate groups in four 5 m 3 circular recycling tanks and 
acclimated over one month at 16°C [42]. Fish were then 
exposed for 21 days to two different temperature treat- 
ments: 16±0.3°C (control groups) and 6.8±0.3°C (cold- 
exposed groups). Because cold-exposed groups stop eating, 
as naturally occurs during winter season, the control 
groups were fasted during the trial, in order to distinguish 
the single effects of low temperature. Warm and cold water 
temperature conditions were maintained by supplying 
tanks with two separate re-circulating systems equipped 
with mechanical and biological filters. Cold water con- 
ditions were achieved by inducing two consecutive 
temperature drops: water was firstly gradually cooled 
from 16°C to 12°C at a rate of 1°C day 1 and then re- 
duced from 12°C to 6.8°C within 24 h. This temperature 
was then kept for 21 days. The experimental design is 
reported in Additional file 5. 

Tissue collection and RNA extraction 

Ten fish per treatment were sampled at different inter- 
vals during acute (0 hour; 0 h; 6 hours, 6 h; 24 hours, 
24 h), and chronic exposure (21 days, 21d) to 6.8°C (see 
Additional file 5). Fish were firstly anaesthetized in a so- 
lution of tricaine methanesulphonate (300 mg/1; MS 222, 
Finquel) in a 30-litre bucket and then euthanized by sev- 
ering their spinal cord. Liver samples (about 30 mg) 
were collected from each specimen, incubated in RNA 
later at 4°C for 24 h and then stored at -20°C. 

Total RNA was extracted from tissue samples using the 
RNAeasy Mini Kit (Qiagen, Hilden, Germany). RNA con- 
centration was determined using a UV-vis spectropho- 
tometer NanoDrop® ND-1000 (NanoDrop Technologies, 
Wilmington, USA). RNA integrity and quality was evalu- 
ated by means of Agilent 2100 Bioanalyzer (Agilent Tech- 
nologies, Palo Alto, CA) according to Ferraresso et al [29]. 
RNA integrity number (RIN) index, that provides a nu- 
merical assessment of RNA integrity, was calculated for 
each sample. This value facilitates the standardization of 
the quality interpretation; only RNAs with RIN >8 were 
used for further microarray experiments. Once single RNA 
samples were extracted and measured, 3 pools (one pool 
consisting of four individuals, two pools_ of three individ- 
uals) per each sampling-time point were prepared. This 
was aimed at getting a sufficient statistical power which is 
reached by a minimum of three biological replicates per 
analysed point and at reducing microarray analysis cost. 

Microarray analysis and data validation 

To perform DNA microarray analyses the microarray 
platform GPL6467 has been used [29]. The annotation 
process of all transcripts represented in microarray was 
performed as described by Ferraresso et al. (2008). To 
update the previous annotation, new blastn have been 
performed against the cDNA database of Danio rerio, 
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Oryzias latipes, Tetraodon Nigrovidis, Takifugu rubripes 
and Gasterosteus aculeatus (cut off e-value of <L0E - 
10). A final number of 8,425 out of 19,734 (42.6%) have 
been annotated in at least one reference database. 

Sample labelling and hybridization were performed ac- 
cording to the Agilent One-Color Microarray-Based Gene 
Expression Analysis protocol for a total of 24 pools (3 bio- 
logical replicates for each conditions and time point). 
Briefly, for each pool 500 ng of total RNA were linearly 
amplified and labeled with Cy3-dCTP. A mixture of 10 
different viral poly-adenilated RNAs (Agilent Spike-In 
Mix) was added to each RNA sample before amplification 
and labeling, to monitor microarray analysis work-flow. 
Labeled cRNA was purified with Qiagen RNAeasy Mini 
Kit, and sample concentration and specific activity (pmol 
Cy3/(ig cRNA) were measured in a NanoDrop® ND-1000 
spectrophotometer. A total of 1,650 ng of labeled cRNA 
was prepared for fragmentation adding 11 \A 10X Blocking 
Agent and 2.2 \A of 25X Fragmentation Buffer, heated at 
60°C for 30 min, and finally diluted by addition with 
55 (il 2X GE Hybridization buffer. A volume of 100 \i\ of 
hybridization solution was then dispensed in the gasket 
slide and assembled to the microarray slide (each slide 
containing four arrays). Slides were incubated for 17 h 
at 65°C in an Agilent Hybridization Oven, subsequently 
removed from the hybridization chamber, quickly sub- 
merged in GE Wash Buffer 1 to disassembly the slides 
and then washed in GE Wash Buffer 1 for approxi- 
mately 1 minute followed by one additional wash in 
pre-warmed (37°C) GE Wash Buffer 2. 

An Agilent G2565BA DNA microarray scanner was 
used to scan arrays at 5 \im resolution, Feature Extrac- 
tion Software 9.S.1 was then used to process and analyse 
array images. The software returns a series of spot qual- 
ity measures in order to evaluate the goodness and the 
reliability of spot intensity estimates. Spike-in control in- 
tensities {Spike-In Viral RNAs) were used to identify the 
best normalization procedure for each dataset. After cyc- 
lic lowess normalization, spike intensities are expected to 
be uniform across the experiments of a given dataset. 
Filtering and normalization analyses were performed 
using R statistical software [52]. Gene expression data 
were deposited in the GEO database under accession 
numbers GSE51442. 

Principal Component Analysis (PCA) was performed 
considering the entire gene expression dataset using 
TMeV software. SAM statistical tests [53] were carried out 
to identify differentially expressed genes between cold and 
control groups (False Discovery Rate (FDR) equal to 5% 
and fold-change (FC) threshold of log 2 FC =1, which cor- 
responds to FC = 2). A more systematic functional inter- 
pretation of differentially transcribed genes was obtained 
through an enrichment analysis using Database for Anno- 
tation, Visualization, and Integrated Discovery (DAVID) 



software [54,55] considering the GO Biological Process 
Database. In order to obtain IDs compatible with DAVID, 
BLASTN (cut off e-value of 1.0 E-5) was used on 
ENSEMBL database to search for significant matches with 
the transcriptomes of the 5 teleost species (see Additional 
file 6). For annotating sea bream transcripts showing sig- 
nificant matches with more than one species, priority was 
assigned first to G. aculeatus (GA), then T. rubripes (TR), 
T. nigroviridis (TN), O. latipes (OL) and D. rerio (DR). 
Retrieved protein sequences from Ensembl genome 
browser (http://www.ensembl.org/biomart/martview), were 
then employed to carry out a Blastp (cut off e-value of 1.0 
E-3) between teleost homologues of sea bream transcripts 
and human Ensembl protein database. Finally, Ensembl 
Human Gene IDs were obtained from the corresponding 
Ensembl protein entries using the BIOMART data mining 
tool (http://www.ensembl.org/biomart/martview/). 

Homologous genes were found for 7,295 transcripts in 
at least one teleost species, while human homologous have 
been identified for 6,835 transcripts. These identifiers were 
used to define a "gene list" and a "background" in the bio- 
informatic tool DAVID, corresponding to differentially 
transcribed genes and to all the transcripts that were rep- 
resented on the array, respectively. Annotation steps per- 
formed for DAVID functional analyses are summarized in 
Additional file 6. The functional annotation was obtained 
for genes differentially expressed between cold and control 
groups setting DAVID for gene count = 2 and ease = 0.1. 
Gene set enrichment analysis (GSEA) was performed on 
the entire gene expression dataset using the Gene Set En- 
richment Analysis v2.0.13 software [56] downloaded from 
the Broad Institute [57]. GSEA analysis was performed by 
using gene symbols retrieved by blastx against UniProt 
database. For enrichment analysis REACTOME gene sets 
were downloaded from the C2 collections in MsigDB v3.1 
(Molecular Signature Database). Pathway analysis was per- 
formed for each gene set independently. Signal2Noise 
metric was employed for gene ranking, and 1000 permuta- 
tions were applied for pathway p-value assignment. 

In order to validate microarray data, ten target genes 
(see Additional file 7) were selected for real-time RT-PCR 
analysis. qPCR assays were performed in a total of 12 sam- 
ples (control and cold-exposed groups at 24 h and 2 Id), 
the same employed for microarray experiments. For each 
selected target gene and for the reference gene (RPL13a, 
Ribosomal Protein L13a), a qRT-PCR assay was designed 
using the SYBRGreen system of detection. Gene-specific 
primers were defined for each transcript with the Primer3 
software. To design intron-spanning primers, putative 
intron-exon boundaries were inferred by comparison with 
homologs of sea bream genes present in high-quality draft 
genome sequences from other fish species {Tetraodon 
nigroviridis, Danio rerio, Oryzias latipes, Takifugu rubripes 
and Gasterosteus aculeatus). 
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One microgram of total RNA for each sample was re- 
verse transcribed to cDNA using Superscript II (Invitro- 
gen™). An aliquot (2.5 |il) of diluted (1:40) cDNA template 
was amplified in a final volume of 10 [A, containing 5 [A 
of mastermix SYBRGreen® Master 2X (Invitrogen) and 
0.25 \A of each gene-specific primer (10 |iM). The ampli 
fication protocol consisted of an initial step of 2 min at 
50°C and 10 min at 95°C, followed by 45 cycles of 10 s at 
95°C and 30 s at 60°C. All experiments were carried out in 
a LightCycler® 480 (Roche Diagnostics). To evaluate the 
efficiency of each assay, standard curves were constructed 
amplifying two-fold serial dilutions of the same cDNA 
which was used as calibrator (Sa_CAL). For each sample, 
the Cp (Crossing point) was used to determine the relative 
amount of target gene; each measurement was made in 
duplicate, and normalized to the reference gene (RPL13a), 
which was also measured in duplicate. RPL13a was chosen 
as reference gene in qRT-PCR assays as it is considered a 
housekeeping gene, and it did not exhibit any significant 
change in microarray data between exposed to low 
temperature and control individuals. To validate pool 
homogeneity and exclude particular differences between 
single samples that constitute the same pool, these were 
individually tested in real-time RT-PCR. In addition, to 
make sure that all targets genes were expressed in the 
calibrator, Sa_CAL, the internal control for each ampli- 
fication reaction, was constructed pooling the same 
amount of three RNA liver samples for each time point 
and condition. 

A non parametric Spearman rank-correlation test was 
applied for each target gene on all 12 samples employed 
for microarray experiments in order to assess the correl- 
ation between the expression values measured respect- 
ively with real-time RT-PCR and microarray probes. 
Statistical analyses were performed using SPSS ver. 12.0. 

Post-mortem examination 

Post-mortem examination was performed on fish experi- 
mentally exposed to low temperature and controls to 
evaluate growth and liver steatosis (fatty liver). Body 
weight (g) and total length (cm) were recorded at the be- 
ginning and end of trial. 

The severity of steatosis was evaluated based on liver 
appearance, colour and consistency [58]. Liver was fur- 
ther weighed and hepatosomatic index (HSI) calculated 
as HSI = LW/BW*100, where LW is the liver weight (g) 
and BW is the body weight. 

The Students t-test was applied to determine the ef- 
fects of cold exposure on growth, LW and HSI and the 
X2-test to test the frequency of liver steatosis 

Data accessibility 

Gene expression analyses were performed using the 
Agilent-016251 Spams aurata oligo microarray (GEO 



accession: GPL6467). Microarray raw and normalized 
fluorescence values were deposited in the GEO database 
(http://www.ncbi.nlm.nih.gov/geo) under accession num- 
ber GSE51442. 



Additional files 



Additional file 1: Differentially expressed genes found between 
control and cold-exposed gilthead seabream at 0 h, 6 h, 24 h and 

21 d. Probes reported in red were over-expressed in cold-exposed 
samples. Probes reported in green were under-expressed in 
cold-exposed samples. Fold change, q-value (%)are also reported. For 
each significant probes has been reported the annotation considering 
amino-acid non redundant database (cut off e-value of <1.0E-10), and 
Ensembl protein database of Donio rerio, Oryzias latipes, Tetroodon 
Nigrovidis, Takifugu rubripes and Gasterosteus aculeatus. 

Additional file 2: Enriched Biological Processes (BP) among 
differentially expressed genes at 0 h, 6 h, 24 h and 21 d. At 24 h 

and 21 d, enrichment analysis was carried separately for over- and 
under-expressed DEGs. 

Additional file 3: Gene set enrichment analysis (GSEA) of DEGs 
found between control and cold-exposed gilthead seabream at 0 h, 
6 h, 24 h and 21 d. 

Additional file 4: Table A. Correlation between microarray and 
real-time RT-PCR expression data, using Spearman rank-correlation test. 
Table B. Comparison of gene expression values between qRT-PCR and 
microarray probes for selected target genes. 

Additional file 5: Experimental design representing fish 
acclimatization and sampling times. 

Additional file 6: Information about the annotation of gilthead 
seabream transcriptome performed in order to obtain IDs 
compatible for the functional analyses. Blastn (cut off e value of <1 OE - 
10) have been performed against the cDNA database of Danio rerio (DR), 
Oryzias latipes (OL), Tetraodon Nigrovidis (TN), Takifugu rubripes (TR) and 
Gasterosteus aculeatus (GA). Retrieved protein sequences from Ensembl 
genome browser, Blastp was then carried out (cut off e-value of 1.0 E-3) 
between teleost homologues of sea bream transcripts and human 
Ensembl protein database. 

Additional file 7: Primer sequences of the ten transcripts analysed 
by real-time RT-PCR to validate Microarray data. 26S: 265 proteasome 
complex subunit; FAS: Fatty acid synthase; GPAT: Glycerol-3-phosfateacyl 
transferase; PKC: Protein kinase C; KAT: 3-ketoacylCoA thiol ase; ILF2: 
lnterleukin-2; MDH: Malate dehydrogenase; SOD: Superoxide dismutase; 
CD59; ACP; RPL13a: Ribosomal Protein LI 3a. 
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HSI: Hepatosomatic index; FC: Fold change; FDR: False Discovery rate; 
DEGs: Differentially expressed genes; HMOX1: Heme oxygenase (decycling) 1; 
CEBPG: CAAT/enhancer binding protein gamma; DDIT3/CHOP: DNA-damage- 
inducible transcript; UPR: Unfolded Protein Response. KEAP1, Kelch-like 
ECH-associated protein 1; NFE2L2/NRF2: Nuclear factor erythroid 2-like 2; 
GPAT-1: Glycerol-3-phosphate acyltransferase mitochondrial precursor; 
ACSL4: Acyl-CoA synthetase long-chain family member 4; 
NAFL: Nonalcoholic fatty liver; ACADS: Acyl-CoA dehydrogenase C-2 to C-3 
short chain; ACADL: Acyl-CoA dehydrogenase long chain; 
HSD17B4: Hydroxysteroid (17-beta) dehydrogenase 4; PPARD: Peroxisome 
proliferator-activated receptor delta; PPARG: Peroxisome proliferator-activated 
receptor gamma; PPARA: Peroxisome proliferator-activated receptor alpha; 
FABP6: Fatty acid binding protein 6; FABP1: Fatty acid binding protein 1; 
FABP2: Fatty acid binding protein 2; FABP7: Fatty acid binding protein 7; 
BP: Biological process; GADD45B: Growth arrest and DNA damage-inducible 
beta; BCL2L14: BCL2-like 14; CIDEC: Cell death-inducing DFFA-like effector c; 
PTEN: Phosphatase and tensin homolog; TP53BP2: Tumor protein p53 
binding protein 2; CASP6: Caspase 6; FLD: Fatty liver disease; XBP1: X-box 
binding protein 1; ATF4: Activating transcription factor 4; IRE1: Inositol- 
requiring enzyme 1; ERAD: Endoplasmic reticulum associated degradation; 



Mininni et al. BMC Genomics 2014, 15:765 
http://www.biomedcentral.eom/l 47 1 -2 1 64/1 5/765 



Page 11 of 12 



ARS: Aminoacyl-tRNA synthetase; EIF: Eukaryotic translation initiation factors; 
ROS: Reactive oxygen species. S0D1, Superoxide dismutase 1; 
S0D2: Superoxide dismutase 2; GSEA: Gene set enrichment analysis; 
PRDX5: Peroxiredoxin 5; GPX4: Glutathione peroxidase 4; PHB1: Prohibitins 1; 
TUBA: a- tubulins; TUBB: (3-tubulins; SCD: A-9 desaturase; HVA: Homeoviscous 
adaptation; GPAT: Glycerol-3-phosfateacyl transferase; 26S: 26S proteasome 
complex subunit; MDH: Malate dehydrogenase; ACP: Acyl-carrier-protein; 
HMGB1: High mobility group box 1. 
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